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m ; 

CN| \ We present a simple novel method for determining the orbital parameters 



of binary pulsars. This method works with any sort of orbital sampling, no 
matter how sparse, provided that information on the period derivatives is 
available with each measurement of the rotational period of the pulsar, and 



> 
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it is applicable to binary systems with nearly circular orbits. We use the 
. technique to precisely estimate the hitherto unknown orbital parameters of 

o ■ 

■ two binary millisecond pulsars in the globular cluster 47 Tucanae, 47 Tuc S and 

Oh' T. The method can also be used more generally to make first-order estimates 

of the orbital parameters of binary systems using a minimal amount of data. 
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1 INTRODUCTION 

The discovery of a pulsar whose period changes significantly usually indicates that it is 
a member of a binary system. It is then important to determine the Keplerian orbital 
parameters of the system in order to investigate the astrophysics of the two stars and to 
obtain a coherent timing solution for the rotation of the pulsar. The latter provides precise 
astrometric information as well as the period derivative which may give information on 
the magnetic flux density and age of the pulsar or, if it is in a globular cluster, on the 
gravitational field of the cluster ( |Anderson 1992| ; [Phinney 1992| ). 



The usual procedure for obtaining orbital parameters involves fitting a Keplerian model 
to a series of period measurements. Such a procedure works well, provided that it is possible 
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to determine the rotational period on several occasions during a single orbit. However, there 
are often circumstances where this is not the case, such as where interstellar scintillation 
permits only sparse positive detections of a pulsar flCamilo et al. 2000|) . 

We have therefore sought a way of employing the hitherto unused information contained 
in the measured orbital accelerations to obtain a full, coherent solution for the Keplerian 
orbital parameters. 

In what follows, we present a simple procedure for estimating the orbital parameters of 
a binary that is completely independent of the distribution of the epochs of the individual 
observations, provided that the period derivatives are known in each observation. In section 
H], we present the equations for the period and acceleration of a pulsar in an eccentric binary 
system as a function of its position in the orbit. We also provide an analysis of the circular 
case and show how the orbital parameters can be obtained from a plot of the acceleration 
of the pulsars against rotational period. In section |3], we demonstrate the application of 
the procedure to a number of pulsars in the globular cluster 47 Tucanae which have known 
orbital parameters, and then to 47 Tuc S and T, neither of which hitherto had an orbital 
solution. In section |], we show how to improve the estimates of the orbital parameters, and 
proceed to refine the orbits of 47 Tuc S and T. 



2 BINARY ORBITS IN THE ACCELERATION/PERIOD PLANE 

Figure [l] represents the orbit of a binary pulsar (point P) around the centre of mass of 

the system (O) projected onto a plane that contains the direction towards the Earth (i.e., 

perpendicular to the plane of the sky, II), and the line of nodes where the orbital plane 

intersects II. The orbit is an ellipse of eccentricity e and projected semi-major axis ap = 

a' P sin i (i is the unknown inclination of the orbit relative to the plane of the sky), r is the 

vector connecting O to P, and r = |r| is the distance of the pulsar to O, ri = \t\\ is the 

projection of this to the line-of-sight. u is the longitude of periastron and / is the angle of 

the pulsar to the periastron measured at O, also called "true anomaly". According to Roy 

(1988), the equation for r\ as a function of / is given by 

<t\ ■ ( . t\ sin(u; + f) 

n{f) =rsm(u + f) =a P — 1 

1 + e cos / 

The time derivative of r/ is the line-of-sight (or radial) velocity, t>/. Using the results in 
Roy (1988), we obtain: 
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Figure 1. Geometrical parameters for an elliptical orbit. 
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where Pb is the orbital period of the binary. The apparent rotational period of the pulsar 
P as a function of / and the intrinsic rotational period Pq is given by 



P(f) =Po[l + 



vi(fY 



-1/2 



Po 1 + 



viifY 



(3) 



c j \ cr J \ c ) 

if the total velocity of the pulsar, v(f) is small compared to c. 

Differentiating equation |2| in time, we obtain the acceleration along the line of sight Ai 
as a function of / 



Mf) - - [ TbJ l 



2 



'\ + e cos f) 2 sin(u; + /) 



(4) 



Plotting Ai(f) as a function of P(f), we obtain a parametric curve that does not depend 
on time, and therefore does not require the solution of Kepler's equation. This curve is 
illustrated in Figure |2] for a binary pulsar of spin period 10 ms and orbital parameters e = 
0.9, u = 140°, Pb = 1 day and x = ap/c = a' P sini/c = Is. 
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Figure 2. A simulated binary system with an orbital period of 1 day, x = Is, e = 0.9 and ui = 140°. The pulsar has an 
intrinsic rotational period of 10 ms. The dashed line is the acceleration as a function of rotational period measured at the 
barycentre of the solar system. This pulsar was "observed" every ten minutes, each observation is represented as a + sign. The 
zones with a smaller concentration of observations correspond to the periastron, where the pulsar spends relatively little time, 
but where accelerations are largest. The cusp is slightly to the right of apastron, which is in the near straight segment where 
the acceleration is close to zero and where the density of observations is greatest. 

All known binaries containing "true" millisecond pulsars, i.e., those with periods below 
20 ms, have rather circular orbits ( Uamilo 1995 ; Lyne 1995| ). For these, we can set e = uo = 



in the equations above, so that equations || and f| reduce to 
2n 

P(f) = P + Pox— cos f = P + Pi cos / (5) 
and 

4vr 2 

A(f) = — —rxc sin / = -A x sin /. (6) 



P 



The track followed by such pulsars in the period/acceleration space is thus an ellipse 
centred on the point (P , 0) and having as horizontal and vertical semi- axes the values Pi 
and A\ respectively, with the pulsar moving in a clockwise direction. 
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Generally, once an ellipse has been found to fit a set of measurements of barycentric 
period and acceleration for a given pulsar, we can easily recover the two relevant orbital 
parameters in a circular orbit: 

P - Pl 2nC (7) 



X 



PA 2 c 



Using these newly determined orbital parameters, we can calculate the angular orbital 
phase for each fcth data point, i.e. for each pair of acceleration and period measured (A)., 
Pk)- 

fc = arctan(-^^^) (9) 

Since the time Tj~ of each observation is also known, we can determine the time of the 
nearest ascending node, or simply T asCj k for each observation: 

^asc,k = T k — P^Pb (10) 
We can use these values to refine the orbital period (section f|). 



3 APPLICATION OF THE PROCEDURE TO PULSARS IN 47 TUCANAE 

In this section, we first test the procedure described above on a number of known binary 
systems in the globular cluster 47 Tucanae and then apply it to two systems with hitherto 
unknown orbits. 

Between August 1997 and August 1999, the 64-m Parkes radio telescope was used to 
make a total of 108 observations of the pulsars in the globular cluster 47 Tucanae (henceforth 
47 Tuc) at a central frequency of 1374 MHz (A = 21cm). About half of these observations 
were of 4.6 hours duration, the remainder were shorter. These data have been used to make 
accurate timing measurements of the pulsars in 47 Tuc, the long integrations being required 
in order to detect the very faint pulsars (Freire et al. 2000). However, it can also be used 
to search for previously undetected pulsars (|Camilo et al. 2000|) using an analysis which 



searches over a range of different accelerations. 

In a first search, targeted at binaries with periods as short as 90 minutes, each of these 
4.6-hour observations was divided into sixteen 17.4-minute segments, and an acceleration 
search was performed on each segment. The trial accelerations ranged from —30 m s -2 to 30 
m s -2 , with a step of 0.3 m s" 2 . This search was remarkably successful, with the discovery 
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of nine new millisecond pulsars, all of which are binary ( |Camilo et al. 2000| ). This increased 
the total number of known pulsars in 47 Tuc to 20. All have periods between 2 and 8 ms, 
and 13 are members of binary systems. 

The reason that new pulsars are still being found in 47 Tuc is that the pulsars scintillate 
strongly at 1374 MHz, with flux density variations over more than one order of magnitude 
in a few hours. This occasionally increases the flux density of a weak pulsar above the 
sensitivity threshold of the Parkes radio telescope, making it detectable for perhaps one or 
two hours in the 17.4-minute sub-integrations. 

If a newly discovered binary system has an orbital period comparable to the length of 
an observation, then it is possible to make a unique determination of the orbital parameters 
of the pulsar using data from a single observation. However, some of the newly-discovered 
binaries (e.g. 47 Tuc S and T) have longer orbital periods. Thus far, it has been impossible to 
determine the orbital parameters of these binaries, because we cannot choose when to observe 
them - that is primarily determined by scintillation, and detections of these pulsars are very 
rare (47 Tuc S was detected on a single observation, and 47 Tuc T on three observations in the 
17.4-minute search). More importantly, the extremely sparse successful detections typically 
occur once every few hundred orbits. For these two pulsars, the number of detections has been 
increased to 4 and 7 respectively, by dividing the observations into 4 70-minute segments 
and conducting acceleration searches on each of these, so increasing the sensitivity by a 
factor of 2 over the search of Camilo et al. (2000). These detections are summarised in Table 
[l], but unfortunately are still much too sparse to determine a unique orbit from the period 
determinations alone. 

The top diagram in Fig. ^| shows a plot of the measured accelerations as a function of 
barycentric period for 47 Tuc W, a 2.35-ms pulsar in a 3.2-hour binary system. The dashed 
ellipse represents the best fit to the data obtained using the implementation described in 
Appendix A. The fit clearly describes the data well. For this particular pulsar, we do not 
observe large negative accelerations which would correspond to times when the pulsar lies 
behind the companion with respect to Earth. This apparent lack of negative accelerations 
is due to the known eclipses in this system, which cover about 40 % of the whole orbit. 

As well as 47 Tuc W, we have repeated the analysis on 47 Tuc J and 47 Tuc U and 
summarise the results from the fitting procedure in Table 0. Given in the table are the fitted 
parameters which can be compared with those obtained form the full, coherent solutions 
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Pulsar 


T 


P 


A 


(S/N) 




1 1V1 J LJ J 


(ms) 


(m s 2 ) 




47 Tuc S 


50741.644 


2.8304388 


-0.78 ±0.06 


16 




50741.693 


2.8304060 


-0.86 ±0.06 


14 




51000.765 


2.8304720 


0.72 ±0.04 


14 




51000.814 


2.8304981 


0.58 ±0.04 


11 




51002.899 


2.8303046 


0.56 ±0.08 


9 




51216.176 


2.8305188 


-0.50 ±0.08 


9 


47 Tuc T 


50740.654 


7.5878889 


-0.92 ± 0.08 


12 




50740.703 


7.5878334 


-0.24 ±0.08 


13 




50746.626 


7.5884100 


1.68 ±0.04 


12 




50746.675 


7.5885810 


1.67 ±0.06 


19 




50981.811 


7.5878836 


0.72 ±0.06 


9 




50983.932 


7.5878552 


-0.44 ± 0.08 


11 




51005.893 


7.5891029 


0.52 ±0.08 


13 




51040.857 


7.5891350 


0.00 ± 0.04 


66 




51335.904 


7.5891358 


0.00 ± 0.02 


11 



Table 1. Measured periods and accelerations for 47 Tuc S and T. S is detected on four different observations of 47 Tuc, T is 
detected on 7 different observations. This represents a total of 6 and 9 measurements of period and acceleration respectively. 
The signal-to-noise ratio (S/N) varies due to scintillation. 



which are of course several orders of magnitude more precise. The fitted values are all 
consistent with these within the quoted errors. 

The other two diagrams in Fig. ||| show the data presented in Table [I]. These few and 
very sparse detections clearly delineate well-defined ellipses and have proved enough to make 
preliminary determinations of the orbits of both of these pulsars, which we made using the 
method outlined in section and Appendix A. The orbital parameters resulting from the 
best fitting ellipses to the data from these two pulsars are given in Table 



4 REFINING THE EPHEMERIS 

As mentioned at the end of section 0, it is possible to refine the orbital period further using 
the local values of the T asCj k, in the same manner that a set of times-of-arrival of pulses from 
a pulsar (TOAs) can be used to improve the estimate of the rotational period of an isolated 
pulsar, once a reasonable estimate of the rotational period is available. This also allows the 
determination of a more precise value for a global T asc . 

As in the pulsar case, the problem is that we don't know precisely how many orbits have 
occurred between any two observations. With data as sparse as that of 47 Tuc S and T, 
this task is normally impossible. However, we already know approximate orbital periods for 
the binary systems from the procedure outlined in section ||] and implemented in Appendix 
A. We also know more than two times of ascending node for each binary, calculated using 



equation 10. Thus, all we have to do is to estimate a binary period and global T asc that 
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Figure 3. Observed accelerations plotted against barycentric spin period for 47 Tuc W (top), 47 Tuc S (middle) and 47 Tuc T 
(bottom), three binary millisecond pulsars in the globular cluster 47 Tuc. The latter two plots were obtained from the numerical 
data in Table |l| and the dashed ellipses represent the best fits. 
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produce the observed T asc k , with the restriction that the orbital period should be consistent 
with that determined in the previous section. 

We have developed two independent methods to perform this task. In the first method, 
we essentially treat each T asc k derived from Eqn. ( [T0| ) as a TOA. For each observation, k, 
we compute a "residual" 

D -^~asc,k -^~asc,l i -^~asc,k -^~asc,l i /i i \ 

R k = p L p J ( U ) 

*B *B 

(where [a\ is the largest integer smaller than a) and plot it against each T asC) k. If the previous 



estimate of Pb used in eq. [Ll] is correct, then all should be randomly distributed around 
zero within the errors. An error in Pb results in a slope in the data points. Hence, we model 
this by performing a least-squares fit of a straight line and repeat this procedure for a set 
of trial Pgs. We search a range of Pbs centred on the previous best estimate, such that the 
total phase difference over the whole time span, T asc n — does not exceed unity (as we 
can assume that our previous guess was sufficiently close to the correct value). We note that 
for each trial Pb, we re-calculate the T asc k 's using the phases, from Eqns. (|9]) and (|T0|). 
The results of the straight line fit are then used to correct the previous estimates for T asc 
and Pb- If the fitted Pb and T asc are at some point coincident with those of the binary, we 
should expect to see a sharp decrease in the residuals given by eq. O 

In the second method, we make direct use of the phase information, 0^, by fitting the 
function 



T — T T — T 

V Pb L Pb J 1 ' 

to our data set of observing epochs and derived phase information (Tk,<pk)- We perform 

a least squares minimization, where due to the non-linear nature of Eqn. [12] we employ a 

robust Downhill-Simplex algorithm ( Press et al. 1992 ) to obtain T asc and Pb- We start the 

simplex routines for various starting points using T asc l and a variation of Pb as in the first 

method. In order to obtain reliable error estimates, we again use a Monte-Carlo simulation 

as described in Appendix A, generating a synthetic data set from the real (T^, </>fc)'s and the 

error estimates for the phase as derived from the first Monte-Carlo application. 

Our experience shows that both methods yield very similar results, usually agreeing 

within the errors. However, depending on the spread of the data points, length of total time 

covered and uncertainties, one or the other method may yield superior results. While the 

orbital solution for 47 Tuc T was refined using the first method, the second method provided 
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Parameter J U W S T 



pa 
1 q 


(ms) 


2.1006335(6) 


4.342831(1) 


2.3523445(8) 


2.830407(8) 


7.588476(4) 


ah 
p o 


(ms) 


2.1006336 


4.342827 


2.3523445 


2.830406 


7.588481 


pa 
1 B 


(d) 


0.119(3) 


0.41(2) 


0.1322(16) 


1.24(14) 


1.12(3) 


pb 
1 B 


(d) 


0.121 


0.43 


0.1330 


1.201724(3) 


1.13 


X a 


00 


0.0396(13) 


0.49(4) 


0.240(3) 


0.79(14) 


1.33(4) 


X b 


(8) 


0.0404 


0.53 


0.243 


0.7659(14) 


1.34 


T 




51000.7979 


51003.0587 


51214.9496 


51000.9649 


51000.3173 


N 




31 


13 


9 


6 


9 



Table 2. Comparison of the orbital parameters for several binaries in the globular cluster 47 Tuc. The parameters indicated 
by an "a" were determined by the method presented in this paper, the parameters indicated by a "b" were determined by the 
timing analysis. The latter are far more accurate (Freire et al. 2000) and their true value is within half of the last significant 
digit quoted. For 47 Tuc S, the "b" parameters are still derived from the period analysis. The method presented in this paper is 
accurate within the error bars derived from Monte Carlo simulations (Appendix A). N is the number of individual measurements 
of acceleration and period; for 47 Tuc J, U and W these were randomly selected from the many measurements of these pulsars 
made by the 17.4-minute search; for 47 Tuc S and T the data used is in Table |l]. 

better results in the case of 47 Tuc S. Applying both methods at the same time serves in 
fact as a good indicator whether the obtained solution is unique. 

With these new, improved Pb and T asc , the orbital parameters can be refined even further, 
using conventional fitting to period and, after that stage, more precise periods derived from 
Times-of Arrival of pulses (see Freire et al. 2000, section 2). In section 3 of that paper we 
present a phase-coherent timing solution for 47 Tuc T; for 47 Tuc S the best solution is: 

P = 2.830406030(3) ms, 

P B = 1.20172465(7) days, 

x = 0.766338(9) s, 

T asc = 51000.964588(8) (MJD). 

No eccentricity has been detected. 

5 CONCLUSION 

The procedure outlined in section |2| and Appendix A is useful even when the binary's orbit 
is well-sampled, because it provides automatically a "first guess" of the orbital parameters, 
which can then be used in conventional period analysis. The quality of this guess obviously 
depends on phase coverage and the precision of the measurement of accelerations. The 
procedure presented in this paper is, however, the only one capable of giving any estimates 
of the orbital parameters when the data samples the orbital period sparsely. For this purpose, 
the epochs of the individual period/acceleration measurements are completely irrelevant as 
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long as the measurements have reasonable orbital phase coverage; their values are only 
needed to make the appropriate period corrections to the barycentre of the solar system. 
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APPENDIX A 

Fitting an ellipse to a set of points can be quite cumbersome, at least from a com- 
putational point of view. However, eqns. (|5]) and (||) can be combined to derive a simple 
expression, which is linear in the unknown parameters and hence easy to fit. By squar- 
ing each of the equations, multiplying one of the squares by an appropriate constant and 
summing them, we obtain the simple equation of a parabola, given in the form 

A 2 = a 2 P 2 + ciiP + a , (1) 

where A and P represent acceleration and period. It is straightforward to model the set of 
squares of the observed accelerations and the corresponding periods, (A 2 ,, P k ) by applying 
a linear least squares fitting algorithm. We use the Singular Value Decomposition method, 
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Figure 1. Distributions of the parameters obtained by fitting fOOOO synthetic data sets generated in the Monte Carlo simulation 
in the case of 47 Tuc W. The values for the orbital parameters so determined are accurate to within 1%. Each distribution 
shows the offset from the best fit value. 

as for instance described by Press et al. (1992), to determine the coefficients ao, a\ and a 2 - 
Q The period and the orbital parameters are then given by 

P>- — (2) 



Pi 



a i 
2nc 



B 



X 



PoV- a 2 

Pb ~~ 



P 2 a 



(3) 



(4) 



2ttPo V fl 2 

Despite the simplicity of the problem, there are large covariances among the fitted pa- 
rameters ao, «i and a 2 , making it difficult to obtain reliable error estimates for the derived 
parameters. For this reason, we employed a Monte-Carlo method in our fitting algorithm, 
which provides not only precise values of the orbital parameters but also reliable and realistic 
estimates of their uncertainties. To achieve this, we generate a total of 10000 synthetic data 
sets, each having the same number of observations as the original set of real observations. 
For a given synthetic data set, the generated random periods and accelerations are Gaussian 
distributed while centred on the real observed data points, with dispersions given by the 
observational uncertainties. 

We perform the described fit algorithm for each synthetic data set and obtain probability 
distributions for the orbital parameters as shown in Figure |l] for the case of 47 Tuc W. We 
choose the maximum of each probability distribution as our best fit value of a given orbital 

* We note as a technical detail, that instead of using the original periods, P^, we use in fact a scaled version of those during the 
fit to avoid numerical problems, i.e. we use an expression C ■ (Pfc — P), where C is a constant and P the mean of all measured 
P k 's. 
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parameter and calculate its uncertainty from the spread in its distribution. The reliability 
of this method is demonstrated in Table where we compare the parameters obtained using 
this method (columns 2 to 4) with the parameters as determined in the subsequent timing 
analysis (columns 5 to 8). The timing method provides far more precise values, however 
the comparison clearly demonstrates that our new method provides accurate values and 
calculates reasonable error estimates. 
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